function wheeler(filenumber)
figure(1); clf; hold on;
timehist=[];


    N=readdata(filenumber,1);
    S=readdata(filenumber,2);
    T=readdata(filenumber,3);
%     multigridlevels=readdata(filenumber,4);
%     bound=readdata(filenumber,5);
%    time=readdata(filenumber,6);
%     dtnom=readdata(filenumber,7);
    dx=readdata(filenumber,8);
    dz=readdata(filenumber,9);
%    mx=readdata(filenumber,10);
%    mz=readdata(filenumber,11);
%     mcs=readdata(filenumber,12);
%     mct=readdata(filenumber,13);
%    mtype=readdata(filenumber,14);
    %mTemp=readdata(filenumber,15);
    %mstressxx=readdata(filenumber,16);
    %mstressxz=readdata(filenumber,17);
    %mstrainxx=readdata(filenumber,18);
    %mstrainxz=readdata(filenumber,19);
    %mfinitestrain=readdata(filenumber,20);
%     vx=readdata(filenumber,21);
%     vz=readdata(filenumber,22);
%     P=readdata(filenumber,23);
Nsurfaces=readdata(filenumber,24);
surfaces=readdata(filenumber,25);

for count=1:filenumber
    timehist(end+1)=readdata(count,6);
    if mod(count,100)==0
        count
    end
end
wheelerdata=zeros(Nsurfaces-5,S);
wheelertime=zeros(1,Nsurfaces-5);


x=zeros(1,S); z=zeros(1,T);
for s=1:S-1 x(s+1)=x(s)+dx(s); end
for t=1:T-1 z(t+1)=z(t)+dz(t); end

for sur=4:Nsurfaces-2
    thickness=surfaces(:,sur)-surfaces(:,sur+1);
    wheelerdata(sur-3,:)=thickness;
    wheelertime(sur-3)=0.5*(timehist(sur-2)+timehist(sur-3));
end

wheelertime=wheelertime/(60*60*24*365*1e6);
[WHEELERTIME X]=meshgrid(wheelertime,x);
clf; pcolor(X,WHEELERTIME,wheelerdata'); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar;
end